Understanding Two-Way Fixed Effects Models

Understanding Two-Way Fixed Effects Models

In my previous post on Fixed vs Between vs Random Effects, I explained three different approaches to panel data analysis. The core strength of fixed effects models is their ability to “control for all time-invariant characteristics of counties.” However, this approach addresses only one dimension of potential confounding in panel data—it misses the crucial time dimension.

Two-way fixed effects (TWFE) models extend one-way fixed effects by controlling for both entity-specific and time-specific unobserved factors. This dual control makes TWFE particularly powerful for causal inference in policy evaluation.

1. What Are Two-Way Fixed Effects?

The Core Idea

While one-way fixed effects control for time-invariant characteristics of counties, two-way fixed effects simultaneously control for:

1. County fixed effects (\(\alpha_i\)): Unobserved county characteristics that don’t change over time
2. Time fixed effects (\(\lambda_t\)): Unobserved factors that affect all counties in the same year

Mathematical Model

\[Y_{it} = \beta X_{it} + \alpha_i + \lambda_t + \epsilon_{it}\]

Where:

• \(Y_{it}\) = outcome for county \(i\) in year \(t\)
• \(X_{it}\) = treatment/explanatory variable
• \(\alpha_i\) = county fixed effects (county-specific intercepts)
• \(\lambda_t\) = time fixed effects (year-specific intercepts)
• \(\epsilon_{it}\) = random error

Why Both Types Matter

County Fixed Effects Control For:

• Geography and climate
• Healthcare infrastructure
• Political culture
• Demographic composition

Time Fixed Effects Control For:

• Economic recessions
• Federal policy changes
• Technological advances
• National health trends

# Load required packages 
library(plm)
library(fixest)
library(dplyr)
library(ggplot2)
library(kableExtra)
library(patchwork)

2. Simulating Data for TWFE

Let’s create a simple county-level dataset with clear county and time effects, using 12 counties observed over 8 years (2015-2022).

Dependent variable: Percentage of uninsured population in each county

Independent variable: Medicaid expansion indicator (binary treatment variable)

• Some counties expand Medicaid in 2018, others never expand
• This creates a staggered treatment design ideal for demonstrating TWFE

Control variable: County unemployment rate

• Helps control for economic conditions that might affect insurance coverage

County fixed effects: Represent persistent differences between counties due to factors like:

• Healthcare infrastructure and provider availability
• Political culture and policy preferences
• Demographics and income levels
• Geographic accessibility to care

Time fixed effects: Capture year-specific shocks affecting all counties equally, such as:

• Economic cycles and recessions
• Federal healthcare policy changes
• National trends in healthcare costs
• Implementation of the Affordable Care Act provisions

This setup allows us to demonstrate how TWFE isolates the causal effect of Medicaid expansion by controlling for both stable county characteristics and common time trends that could otherwise confound our estimates.

set.seed(123)
n_counties <- 12
n_years <- 8
total_obs <- n_counties * n_years

# Generate identifiers
county_id <- rep(1:n_counties, each = n_years)
year <- rep(2015:2022, times = n_counties)

# Simple county fixed effects (different baseline levels)
county_effects <- rep(rnorm(n_counties, mean = 0, sd = 20), each = n_years)

# Simple time fixed effects (common trends affecting all counties)
time_effects_values <- c(-5, -3, 0, 2, 3, 1, -2, 4)  # One for each year
time_effects <- rep(time_effects_values, times = n_counties)

# Treatment: Medicaid expansion (some counties expand in 2018, others never)
treatment_year <- rep(sample(c(2018, NA), n_counties, 
                           replace = TRUE, prob = c(0.5, 0.5)), each = n_years)
medicaid_expansion <- ifelse(is.na(treatment_year), 0, 
                           ifelse(year >= treatment_year, 1, 0))

# Control variable
unemployment <- rep(rnorm(n_counties, mean = 5, sd = 1), each = n_years) + 
               rnorm(total_obs, 0, 0.5)

# Outcome: Uninsured rate
treatment_effect <- -4  # True effect: 4 percentage point reduction
uninsured_rate <- 12 +  # Baseline
                 treatment_effect * medicaid_expansion +  # Treatment effect
                 0.3 * unemployment +  # Unemployment effect
                 county_effects +  # County differences
                 time_effects +   # Time trends
                 rnorm(total_obs, 0, 1)  # Noise

# Create dataset
twfe_data <- data.frame(
  county_id = factor(county_id),
  year = year,
  medicaid_expansion = medicaid_expansion,
  uninsured_rate = uninsured_rate,
  unemployment = unemployment,
  treatment_year = treatment_year
)

# Show sample
head(twfe_data, 12) %>%
  kable(digits = 2, caption = "Sample Two-Way Fixed Effects Data") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Sample Two-Way Fixed Effects Data

county_id	year	medicaid_expansion	uninsured_rate	unemployment	treatment_year
1	2015	0	-2.20	5.91	2018
1	2016	0	-0.47	5.55	2018
1	2017	0	1.67	6.15	2018
1	2018	1	0.56	6.14	2018
1	2019	1	3.07	6.11	2018
1	2020	1	0.06	6.05	2018
1	2021	1	-3.37	5.98	2018
1	2022	1	2.07	5.67	2018
2	2015	0	1.66	4.37	2018
2	2016	0	6.83	4.34	2018
2	2017	0	7.19	4.18	2018
2	2018	1	7.46	4.42	2018

3. Visualizing the Need for TWFE

The visualization shows three key patterns:

# Add treatment group labels
twfe_data <- twfe_data %>%
  mutate(treatment_group = ifelse(is.na(treatment_year), "Never Treated", "Treated 2018"))

# Plot 1: Raw trends show both county differences and time patterns
ggplot(twfe_data, aes(x = year, y = uninsured_rate, color = treatment_group)) +
  geom_line(aes(group = county_id), alpha = 0.5) +
  stat_summary(fun = mean, geom = "line", size = 2) +
  geom_vline(xintercept = 2018, linetype = "dashed", color = "black") +
  scale_color_manual(values = c("Never Treated" = "blue3", "Treated 2018" = "red3")) +
  labs(title = "Raw Data: County Trends Over Time",
       subtitle = "Shows county differences, time trends, and treatment effects mixed together",
       x = "Year", y = "Uninsured Rate (%)", color = "Group") +
  theme_minimal()

1. 1st Figure: Individual county trends vary, with treatment and control groups following different paths

# Plot 2: Average by year (shows time effects)
time_means <- twfe_data %>%
  group_by(year) %>%
  summarize(avg_uninsured = mean(uninsured_rate))

ggplot(time_means, aes(x = year, y = avg_uninsured)) +
  geom_line(size = 1.5, color = "darkgreen") +
  geom_point(size = 3, color = "darkgreen") +
  labs(title = "Average Uninsured Rate by Year",
       subtitle = "Shows common time trends affecting all counties",
       x = "Year", y = "Average Uninsured Rate (%)") +
  theme_minimal()

2. 2nd Figure: Common time patterns affect all counties (economic cycles, federal policies)

# Plot 3: Average by county (shows county effects)
county_means <- twfe_data %>%
  group_by(county_id) %>%
  summarize(avg_uninsured = mean(uninsured_rate)) %>%
  mutate(county_rank = rank(avg_uninsured))

ggplot(county_means, aes(x = county_rank, y = avg_uninsured)) +
  geom_col(fill = "steelblue", alpha = 0.7) +
  labs(title = "Average Uninsured Rate by County",
       subtitle = "Shows persistent differences between counties",
       x = "County (Ranked by Uninsured Rate)", y = "Average Uninsured Rate (%)") +
  theme_minimal()

3. 3rd Figure: Persistent differences between counties (some always higher/lower than others)

4. The TWFE Demeaning Process: Step by Step

The key to understanding TWFE is seeing how it removes both county and time effects through “demeaning”:

# Step 1: Calculate means
twfe_data <- twfe_data %>%
  group_by(county_id) %>%
  mutate(county_mean = mean(uninsured_rate)) %>%
  ungroup() %>%
  group_by(year) %>%
  mutate(year_mean = mean(uninsured_rate)) %>%
  ungroup() %>%
  mutate(
    overall_mean = mean(uninsured_rate),
    # Step 2: Apply two-way demeaning transformation
    uninsured_demeaned = uninsured_rate - county_mean - year_mean + overall_mean
  )

The Mathematics of Two-Way Demeaning

The two-way fixed effects transformation follows this formula:

\(Y_{it}^* = Y_{it} - \bar{Y}_i - \bar{Y}_t + \bar{Y}\)

Where:

• \(Y_{it}^*\) = the demeaned (transformed) outcome
• \(Y_{it}\) = the original outcome for county \(i\) in year \(t\)
• \(\bar{Y}_i\) = the average outcome for county \(i\) across all years
• \(\bar{Y}_t\) = the average outcome in year \(t\) across all counties
  
• \(\bar{Y}\) = the overall average across all counties and years

Understanding the Logic Step by Step

Step 1: Why subtract the county mean (\(\bar{Y}_i\))?

# Visualization of each step
ggplot(twfe_data, aes(x = year, y = uninsured_rate, color = treatment_group)) +
  geom_line(aes(group = county_id), alpha = 0.6) +
  stat_summary(fun = mean, geom = "line", size = 2) +
  scale_color_manual(values = c("Never Treated" = "blue3", "Treated 2018" = "red3")) +
  labs(title = "Step 1: Original Data",
       subtitle = "Raw uninsured rates with all variation",
       x = "Year", y = "Uninsured Rate (%)", color = "Group") +
  theme_minimal()

Each county has its own baseline level due to permanent characteristics (geography, demographics, institutions). By subtracting each county’s average, we remove these time-invariant differences. After this step, all counties are centered around zero, eliminating the county fixed effects.

Step 2: Why subtract the year mean (\(\bar{Y}_t\))?

# After removing county effects
twfe_data <- twfe_data %>%
  mutate(after_county_demean = uninsured_rate - county_mean + overall_mean)

ggplot(twfe_data, aes(x = year, y = after_county_demean, color = treatment_group)) +
  geom_line(aes(group = county_id), alpha = 0.6) +
  stat_summary(fun = mean, geom = "line", size = 2) +
  scale_color_manual(values = c("Never Treated" = "blue3", "Treated 2018" = "red3")) +
  labs(title = "Step 2: After Removing County Effects",
       subtitle = "County differences eliminated, time patterns remain",
       x = "Year", y = "County-Demeaned Rate", color = "Group") +
  theme_minimal()

Now we can see that the distance between county lines has narrowed considerably! After removing county-specific effects, the data reveals a clear common time pattern across all counties.

This pattern reflects year-specific factors that affect all counties equally - such as economic cycles, federal healthcare policies, national insurance market trends, and macro-level health crises. By subtracting each year’s average from the county-demeaned data, we remove these time-specific shocks, eliminating the systematic year-to-year variation that would otherwise confound our treatment effect estimates.

What remains after this second demeaning step is variation that differs across counties within the same time period - this is the key identifying variation that allows TWFE to isolate causal effects while controlling for both persistent county differences and common temporal trends.

Step 3: Why add back the overall mean (\(\bar{Y}\))?

# After removing both county and time effects
ggplot(twfe_data, aes(x = year, y = uninsured_demeaned, color = treatment_group)) +
  geom_line(aes(group = county_id), alpha = 0.6) +
  stat_summary(fun = mean, geom = "line", size = 2) +
  scale_color_manual(values = c("Never Treated" = "blue3", "Treated 2018" = "red3")) +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.7) +
  labs(title = "Step 3: After Removing County & Time Effects",
       subtitle = "Only treatment effects and idiosyncratic variation remain",
       x = "Year", y = "Two-Way Demeaned Rate", color = "Group") +
  theme_minimal()

Now we see a clearer pattern emerge between never-treated counties and treated counties! The treatment effect becomes much more apparent after removing both county-specific baselines and common time trends.

This transformation works through careful mathematical adjustment: when we subtract both county means (\(\bar{Y}_i\)) and year means (\(\bar{Y}_t\)), we inadvertently subtract the overall mean (\(\bar{Y}\)) twice. This double-subtraction occurs because the overall mean is embedded within both county averages and year averages. Adding back \(\bar{Y}\) corrects this mathematical artifact, ensuring our transformation removes only the county-specific and time-specific components while preserving the interpretable scale of the original data.

The result is data that isolates the treatment effect by comparing how treated counties differ from control counties at the same point in time, after accounting for both persistent county characteristics and common temporal shocks.

What Remains After Demeaning?

After this transformation, only two sources of variation remain:

1. Treatment effects: Changes in treated counties that differ from changes in control counties at the same time
2. Idiosyncratic shocks: Random county-year specific factors (the error term)

This is the “identifying variation” - variation that occurs within counties and within time periods simultaneously. TWFE uses this remaining variation to estimate treatment effects while controlling for all county-specific and time-specific confounders.

Intuitive Example

Consider county A with high baseline uninsured rates and county B with low baseline rates. In a recession year, both counties see increases in uninsured rates.

• County demeaning removes the baseline differences between A and B
• Time demeaning removes the common recession effect
  
• What’s left is only the variation that differs between counties within that same year - which could be due to treatment effects

This is why TWFE provides such powerful control for confounding factors.

5. Estimating TWFE Models

5.1 Within Transformation (Demeaning)

Again, two-way fixed effects estimation follows a sequential demeaning process that differs from conventional one-way fixed effects:

Step 1: Individual demeaning (Same as conventional FE) \[\tilde{Y}_{it} = Y_{it} - \bar{Y}_i\] \[\tilde{X}_{it} = X_{it} - \bar{X}_i\]

This removes county-specific means from both the outcome and explanatory variables.

In our example: We subtract each county’s average uninsured rate and average Medicaid expansion status from their respective yearly values. This eliminates persistent differences between counties - some counties that are always high-uninsured vs. low-uninsured due to demographics, healthcare infrastructure, or political culture.

Step 2: Time demeaning of demeaned variables (Major difference from conventional FE) \[\ddot{Y}_{it} = \tilde{Y}_{it} - \bar{\tilde{Y}}_t\] \[\ddot{X}_{it} = \tilde{X}_{it} - \bar{\tilde{X}}_t\]

This removes time-specific means from the already county-demeaned variables. Note that \(\bar{\tilde{Y}}_t\) is the average of the county-demeaned outcomes in year \(t\).

In our example: We subtract the yearly averages of the county-demeaned data. This removes common time trends affecting all counties - such as economic cycles, federal policy changes, or national healthcare trends that cause all counties’ uninsured rates to move together over time.

Step 3: OLS on double-demeaned variables \[\ddot{Y}_{it} = \beta \ddot{X}_{it} + \ddot{\epsilon}_{it}\]

Run standard OLS regression on the double-demeaned variables. The coefficient \(\beta\) gives us the two-way fixed effects estimate.

In our example: We estimate how changes in Medicaid expansion (relative to county and time averages) affect changes in uninsured rates (relative to county and time averages). This isolates the pure treatment effect by comparing treated counties to control counties at the same point in time, after removing both county-specific baselines and time-specific shocks.

Why This Sequential Process Matters

The key insight is that we don’t just subtract raw year means - we subtract the means of the already county-demeaned data. This ensures we’re removing time effects from data that has already had county effects removed.

In our context: This means we’re not just controlling for the fact that 2020 was a bad year for healthcare (raw time effect), but rather how 2020 affected counties after accounting for their baseline differences. This prevents time effects from being confounded with county-level characteristics.

Mathematically, this is equivalent to the single-step formula we saw earlier: \[Y_{it}^* = Y_{it} - \bar{Y}_i - \bar{Y}_t + \bar{Y}\]

But the two-step process makes it clearer how TWFE builds on conventional fixed effects methodology.

5.2 Model Comparison

Now let’s estimate all four specifications using the fixest package, which provides fast and modern fixed effects estimation:

# Model 1: No fixed effects (pooled OLS)
model_none <- feols(uninsured_rate ~ medicaid_expansion + unemployment, 
                   data = twfe_data)

# Model 2: County fixed effects only
model_county <- feols(uninsured_rate ~ medicaid_expansion + unemployment | county_id,
                     data = twfe_data)

# Model 3: Time fixed effects only  
model_time <- feols(uninsured_rate ~ medicaid_expansion + unemployment | year,
                   data = twfe_data)

# Model 4: Two-way fixed effects (county + time)
model_twfe <- feols(uninsured_rate ~ medicaid_expansion + unemployment | county_id + year,
                   data = twfe_data)

# Display all models in a comparison table
etable(model_none, model_county, model_time, model_twfe,
       headers = c("No FE", "County FE", "Time FE", "Two-Way FE"),
       digits = 3,
       se.below = TRUE)

                       model_none   model_county     model_time     model_twfe
                            No FE      County FE        Time FE     Two-Way FE
Dependent Var.:    uninsured_rate uninsured_rate uninsured_rate uninsured_rate
                                                                              
Constant                  44.4***                                             
                          (8.69)                                              
medicaid_expansion       -11.5**          0.712*      -23.0***        -3.52***
                          (3.53)         (0.273)       (0.727)        (0.334) 
unemployment              -5.13**         0.104        -5.19***       -0.176  
                          (1.82)         (0.736)       (0.846)        (0.166) 
Fixed-Effects:     -------------- -------------- -------------- --------------
county_id                      No            Yes             No            Yes
year                           No             No            Yes            Yes
__________________ ______________ ______________ ______________ ______________
S.E. type                     IID  by: county_id       by: year  by: county_id
Observations                   96             96             96             96
R2                        0.17989        0.98135        0.31740        0.99753
Within R2                      --        0.01276        0.30731        0.42858
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Understanding fixest Syntax

The fixest package uses a clean syntax where fixed effects are specified after the | symbol:

• y ~ x | fe1: One-way fixed effects for fe1
• y ~ x | fe1 + fe2: Two-way fixed effects for fe1 and fe2
  
• y ~ x: No fixed effects (standard OLS)

# Extract and compare treatment effects specifically
models_list <- list(
  "No FE" = model_none,
  "County FE" = model_county,
  "Time FE" = model_time, 
  "Two-Way FE" = model_twfe
)

# Create coefficient comparison dataframe
coef_comparison <- data.frame(
  Model = names(models_list),
  Coefficient = sapply(models_list, function(x) coef(x)["medicaid_expansion"]),
  SE = sapply(models_list, function(x) se(x)["medicaid_expansion"])
) %>%
  mutate(
    CI_Lower = Coefficient - 1.96 * SE,
    CI_Upper = Coefficient + 1.96 * SE,
    True_Effect = -4,  # Our known true effect
    Bias = round(Coefficient - True_Effect, 3)
  )

# Display detailed comparison
coef_comparison %>%
  mutate(across(c(Coefficient, SE, CI_Lower, CI_Upper), ~round(.x, 3))) %>%
  kable(caption = "Treatment Effect Estimates: Medicaid Expansion on Uninsured Rate") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Treatment Effect Estimates: Medicaid Expansion on Uninsured Rate

	Model	Coefficient	SE	CI_Lower	CI_Upper	True_Effect	Bias
No FE.medicaid_expansion	No FE	-11.452	3.529	-18.368	-4.535	-4	-7.452
County FE.medicaid_expansion	County FE	0.712	0.273	0.176	1.248	-4	4.712
Time FE.medicaid_expansion	Time FE	-22.995	0.727	-24.419	-21.570	-4	-18.995
Two-Way FE.medicaid_expansion	Two-Way FE	-3.523	0.334	-4.177	-2.869	-4	0.477

# Visualize the comparison with specified order
coef_comparison$Model <- factor(coef_comparison$Model, 
                               levels = c("No FE", "County FE", "Time FE", "Two-Way FE"))

ggplot(coef_comparison, aes(x = Model, y = Coefficient)) +
  geom_hline(yintercept = -4, color = "grey40", linetype = "dashed", size = 1.2) +
  geom_point(size = 4, color = "red", shape = 15) +
  geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper), 
                width = 0.2, size = 1, color = "black") +
  coord_flip() +
  scale_x_discrete(limits = rev(levels(coef_comparison$Model))) +  # Reverse for coord_flip
  labs(title = "Treatment Effect Estimates Across Specifications",
       subtitle = "Grey dashed line shows true simulated effect (-4 percentage points)",
       x = "Model Specification", 
       y = "Estimated Effect on Uninsured Rate (pp)",
       caption = "Error bars show 95% confidence intervals") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14))

This visualization shows how different model specifications dramatically affect our estimates:

1. No Fixed Effects: Severely overestimates the effect (-12 pp) due to uncontrolled selection bias - counties choosing Medicaid expansion differ systematically from non-expanding counties.

2. County FE Only: Underestimates (-2 pp) by ignoring time trends. Fails to account for broader healthcare improvements happening during the study period.

3. Time FE Only: Dramatically overestimates (-22 pp) by ignoring county differences. Counties with better healthcare systems were more likely to expand, confounding results.

4. Two-Way FE: Provides the most credible estimate (-4 pp), closest to our true simulated effect, by controlling for both county characteristics and time trends simultaneously.

Visualize Fixed Effects

The fixest package allows us to easily extract and examine the estimated fixed effects:

Interpreting the Fixed Effects

county_fe <- fixef(model_twfe)$county_id
# Display county fixed effects
county_fe_df <- data.frame(
  County_ID = names(county_fe),
  County_FE = round(county_fe, 2)
) %>%
  arrange(County_FE)
ggplot(county_fe_df, aes(x = reorder(County_ID, County_FE), y = County_FE)) +
  geom_col(fill = "steelblue", alpha = 0.7) +
  coord_flip() +
  labs(title = "County Fixed Effects",
       subtitle = "Persistent differences between counties (relative to reference)",
       x = "County ID", y = "Fixed Effect Estimate") +
  theme_minimal()

County Fixed Effects: These represent persistent differences between counties after controlling for time trends and other covariates.

• Positive values: Counties with systematically higher uninsured rates
• Negative values: Counties with systematically lower uninsured rates
• These differences capture unobserved county characteristics like healthcare infrastructure, political culture, or geographic factors

time_fe <- fixef(model_twfe)$year

# Display time fixed effects  
time_fe_df <- data.frame(
  Year = names(time_fe),
  Time_FE = round(time_fe, 2)
)

ggplot(time_fe_df, aes(x = Year, y = Time_FE)) +
  geom_col(fill = "darkgreen", alpha = 0.7) +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.7) +
  labs(title = "Time Fixed Effects", 
       subtitle = "Year-specific shocks affecting all counties",
       x = "Year", y = "Fixed Effect Estimate") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Time Fixed Effects: These capture year-specific shocks that affected all counties equally.

• They represent common trends like economic cycles, federal policy changes, or national health initiatives
• The pattern shows how macro-level factors influenced uninsured rates across all counties in each year

Key Insight: By estimating and removing these fixed effects, our treatment effect estimate is based only on variation that differs across counties within the same time period - providing cleaner causal identification.

Key Insights

The results show how different specifications can lead to different conclusions:

1. No Fixed Effects: Likely biased due to omitted county and time factors
2. County FE Only: Controls for county differences but not time trends
3. Time FE Only: Controls for time trends but not county differences
   
4. Two-Way FE: Most credible estimate, controlling for both sources of bias

The two-way fixed effects estimate should be closest to our true simulated effect of -4 percentage points.